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Abstract 

A particle in the Henon-Heiles potential can escape when its energy is above the threshold value 
Eth = g- We report a theoretical study on the the escape rates near threshold. We derived an 
analytic formula for the escape rate as a function of energy by exploring the property of chaos. 
We also simulated the escaping process by following the motions of a large number of particles. 
Two algorithms are employed to solve the equations of motion. One is the Runge-Kutta-Fehlberg 
method, and another is a recently proposed fourth order symplectic method. Our simulations show 
the escape of Henon-Heiles system follows exponential laws. We extracted the escape rates from 
the time dependence of particle numbers in the Henon-Heiles potential. The extracted escape rates 
agree with the analytic result. 
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I. INTRODUCTION AND MODEL 

Bauer and BertschjjJ studied the decay laws of chaotic and non-chaotic billiards with 
windows. They found the number of classical particles remaining inside chaotic billiards 
decreases exponentially, but no-chaotic billiards decay according to power laws. The results 
suggest that the exponential law is connected to the chaotic dynamics. However, an excep- 
tion is the circular biUiard, which is iutegrable but decays expoueutiaLyfl. Expetirueuta. 
studies on the decay laws of an elbow cavity using microwaves have also been reportedp(. 
We consider the Henon-Heiles systemj^j, with the following Hamiltonian 

H = ~(pl+pl) + U(x,y), (1) 
U(x,y) = ^(x 2 + y 2 ) +x 2 y -^y 3 , 

where x and y are the coordinates, p x and p y are the momenta. The mass of the particle 
is set to one for convenience. This system exhibits both regular motion and chaotic motion 
depending on the energy of the system, and it has been studied from statistical, semiclassical 
and other perspectives E, s[ O, Q. Recently Brack et al have calculated the density of 
states above threshold [111], but the escape problem has not been addressed so far. 

Numerical studies show the motion of Henon-Heiles system is regular for E < 1/12. 
When E is greater than 1/12, the fraction of chaotic region in phase space increases with 
increasing energy until E = 1/6 the whole phase space is chaotic. E t h = 1/6 is the threshold 
energy of this system. When E > E t h, a particle in the potential well can escape. Fig. 1 
shows the contours of the potential U(x, y). There are three saddle points P\{x — 0, y — 1), 
P 2 (-VS/2, -1/2) and P 3 {V3/2, -1/2). All contours with energy less than 1/6 are closed. 
A particle with energy less than 1/6 always moves inside the closed contour and it remains 
in the well. The contour with E = 1/6 is the equilateral triangle P1P2P3. The contours with 
energies larger than 1/6 are not closed. There are three openings at the three saddle points. 
A particle with energy above 1/6 can escape from the well via the three openings. Because 
this system is chaotic above threshold, the escape in the Henon-Heiles system should also 
follow an exponential law. Assuming N(0) random particles with the same energy in the 
Henon-Heiles well at t = 0, the number of particles at t will be 

N(t) = N(0)exp(-oct), (2) 
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where a is an energy dependent decay rate. The purpose of this article is to verify Eq.(2) 
and to estimate a for different energies. 



II. AN ESCAPE RATE FORMULA 

We can derive a formula for the escape rate as a function of energy above threshold by 
using chaotic property of the Henon-Heiles system. We draw a line perpendicular to the 
escape direction through each saddle point, they are line A\B\,A 2 B2 and A 3 B 3 in Fig. 1. 
For any energy E above threshold, we define the potential well of the Henon-Heiles system to 
be the region restricted by the three disconnected contour lines and the three straight lines 
A\B\,A2B2 and A3B3. The motion in the well is assumed to be ergodic because of chaos. 
The distribution on the energy shell is given generally by ip(q,p) = j dqdpi{E-H(q pj) ' wnere 
q,p are the coordinates and momenta^]. For our two-dimensional system, it is easy to work 
out the results: the distribution in (x,y) is uniform inside the well, and once the particle's 
position is given, the magnitude of the momentum is fixed by the Hamiltonian in Eq.(l) 
and the direction of the momentum is uniformly distributed in [0, 2tt]. We define the energy 
above threshold AE = E — E^. We use 9 to represent the direction of momentum relative 
to the y axis. We use S(AE) to denote the area of the well. Then the distribution in the 
variables (x, y, 9) can be expressed as p(x, y, 9) = 2n s(AE) • Given N particles in the well, the 
number of particles leaving the well through the opening at the saddle point Pi in unit time 
can be written as N f dx §_^i 2 d9p(x, y, 9)\v(x, y)\cos(9), where the integral in x is along the 
line A1B1 and is restricted to the classical allowed part. We note the three openings of the 
system are symmetric. Therefore the number of particles leaving the well in unit time from 
three openings are just three times of the above result. The change of N with respect to t is 



dN(t) r /2 rV 2AE / 3 , 

—±± = -3N(t)p / cos(9)d9 / ^/2{AE - 3x 2 /2)dx (3) 

dt J-n/2 J-^/2AE/3 

= -2nV3AEpN(t), (4) 

which gives the escape rate a(AE) = . 

There is no analytical formula for the area of the well S(AE). We have applied Monte 
Carlo method to calculate the function of AE. The numerical results are repre- 

sented by dots in Fig.2. We found the numerical results can be represented vary well by the 
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quadratic polynomial 

S{AE) = S + S 1 AE + S 2 {AE) 2 (5) 

where So = is the area of equilateral triangle PiP 2 P 3 . By fitting Eq.(5) to the numerical 
results in Fig. 2 using least squares, we have determined the values for the other two coef- 
ficients, Si = 9.656 and S 2 = —22.61. The line in Fig.2 is the fitted quadratic polynomial. 
We finally have the formula for the escape rate in the Henon-Heiles system, 

= S + Sl AE + S 2 (AEr (6) 
Very close to the threshold, a power expansion of Eq.(6) may be useful. Define 

00 

a{AE) = Y,BiAE\ (7) 
1 

The coefficients Bi can be expressed in S , Si and They are B\ — B 2 — — ^p- 1 , and 

SB I s B 

others can be obtained from the following iteration formula, Bi = — 1 l ~ 1 s 2 l ~' 2 , i = 3, 4, .... 
The numerical values for the first four coefficients in the power expansion of Eq.(6) are 
Bi = 4/3,^2 = -9.9115,53 = 96.8743, and B 4 = -892.547. The first term in the power 
expansion is proportional to AE, it gives the scaling of the escape rate in the Henon-Heiles 
system just above threshold. 



III. NUMERICAL SIMULATIONS 



In our numerical simulations of the escape processes, we follow a large number of particles. 
We monitor the number of particles N(t) remaining in the potential well as a function 
of time and then extract the escape rate. For any energy above threshold, we initially 
place iV(0) particles in the well according to the distribution p(x,y,9) = 2 TrS(AE) • ^ ms 
distribution sets the initial conditions for the particles. The trajectory of each particle is 
then followed by numerically solving the Hamilton's equations. We used two algorithms to 
integrate Hamilton's equation. Runge-Kutta-Fehlberg (RKF) is the first algorithm jljj]. In 
this algorithm the error in each step can be controlled by setting the relative tolerance and 
the absolute tolerance. In all our calculations, we set the absolute tolerance to 10~ 9 . The 
second algorithm (CC) was proposed recently by Chin and Chen[l^], it is a fourth order 
forward symplectic algorithm. It is generally believed that symplectic algorithms are better 
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and can follow the true dynamics longer because they preserve the symplectic structures 
of Hamilton's equation. The explicit algorithm for advancing the system forward from t to 
t + e is 

Pl = p(<) + ieF(q(i)) (8) 

/ ,s 1 

qi = qW + 2 e Pi 

4 ~ . 
p 2 = Pi + gCF(qi) 

q(i + 1) = qi + iep 2 

p(i + l) = P2 + ieF(q(2 + l)). 

Note F = — VC/ and F = F + ^e 2 V(|F| 2 ) includes an correction to the original force. The 
time step size e can be varied to control integration errors. 

In Fig. 3 we show the two trajectories calculated using the two algorithms. The two 
trajectories start from the same initial condition: the position is at (x — 0, y — 0.16), the 
energy is E — 0.18, and the direction of momentum is in the positive x-axis. Fig. 3(a) is 
the trajectory obtained using RKF algorithm with relative tolerance 10~ 8 . Fig. 3(b) is the 
trajectory obtained using CC method with a time step size e = 0.04. We have verified the 
accuracy of both calculations. When the time step size and relative tolerance were reduced 
further, the trajectories did not show noticeable change. Fig. 3(a) and Fig.3(b) show clearly 
the two trajectories obtained using the two algorithms stay close for some time and then 
separate. In Fig. 3(a), the particle escapes at t — 299, while in Fig.3(b), the particle escapes 
at t — 172. The increased separation of the two trajectories calculated with two different 
algorithms starting with the same initial condition reflects the difficulty to follow chaotic 
motions for a long time. Nevertheless we can still extract accurate escape rates as shown 
below. 

A large number of particles are used in the simulations of escape process. For example, 
we initially put N(0) = 15326 particles with the energy AE = 0.0234 above threshold 
according to p(x, y, 6) = 27T s(ae) ^ n ^ e well. We advanced this system by following the 
trajectory of each particle using RKF method with relative tolerance 10~ 8 . We recorded the 
number of particles N(U) remaining in the well in time step At = 0.628. Fig. 4(a) shows 
N(t) as a function of time in log-linear scale. Fig. 4(b), Fig. 4(c) and Fig. 4(d) show similar 
results for different energies. One can see the curves in Fig. 4 in all the cases are almost 
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straight lines, indicating the escape of Henon-Heiles system follows exponential laws. For 
each energy, we can extract the escape rate from the simulated N(t). We used the simulated 
N(t) from time t = to a time when ten percent of the particles have escaped and fitted it 
to lnN(t) = c — at using least squares. The fitted parameter a is the extracted escape rate 
at the corresponding energy. 

We compare in Fig.5 the extracted escape rates and the analytic result in Eq.(6) as a 
function of energy above threshold. The solid line is from the formula in Eq.(6). The circles 
are the numerical results from RKF method with a relative tolerance 10~ 8 , the diamonds are 
the numerical results obtained from CC method (Ref.13) with time step 0.04. We verified 
the numerical accuracy with smaller tolerance and smaller step size, and we found the results 
did not change. The dashed line also shows the result of the power expansion of Eq.(6) with 
only the first four terms. The extracted results from the two algorithms agree well. They 
also agree with the formula in Eq.(6). The first four terms power expansion is accurate close 
to threshold, when AE is greater than 0.05, it starts to deviate from the analytic result in 
Eq(6) and from the numerical results. 

IV. CONCLUSIONS 

We have derived a formula in Eq.(6) for the escape rate of the Henon-Heiles system. We 
also simulated the escape process by following a large number of trajectories. We used a 
symplectic algorithm and a non-symplectic algorithm to advance each particle's trajectory 
in time. We monitored the number of particles remaining in the potential well in time. 
We found the escape follow exponential laws similar to the chaotic billiard systems[lj. The 
extracted escape rates using the two algorithms agree with each other, and and they also 
agree with the analytic rate formula for Henon-Heiles system in Eq.(6). 
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FIG. 1: Equipotential lines of the function U(x,y) in Eq.(l). The points P\,Pi and P3 are saddle 
points. 
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FIG. 2: The area of the well (region bounded by the three disconnected contour lines and the three 
straight lines A\B\,AiB2 and ^3-63) as a function of energy above threshold AE. The dots are 
results calculated using Monte Carlo method. The solid line is the fitted quadratic polynomial in 
Eq.(5). 
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FIG. 3: Trajectories calculated from the two algorithms in the Henon-Heiles system. Both tra- 
jectories start from the position (x = 0, y = 0.16) with energy E = 0.18, the initial direction of 
momentum points to positive x-axis. (a) RKF method (Ref.12) with relative tolerance 10~ 8 and 
time range (0-299). (b) CC method of Chin and Chen (Ref.13) with step 0.04 and time range 
(0-172). 
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FIG. 4: The number of particles in the well as a function of time for four energies above threshold, 
(a) AE = 0.0234, N(t = 0) = 15326;(b)A£ = 0.0534, N(t = 0) = 17392; (c)AE = 0.0734, 
N(t = 0) = 18885;(d)A£' = 0.0934, N(t = 0) = 19942. The dots are numerical results obtained 
using RKF method with a relative tolerance 10~ 8 . 
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FIG. 5: The escape rate a vs AE. The solid line is the formula in Eq.(6). The circles are the 
numerical results from RKF method with a relative tolerance 10 -8 , and the diamonds are the 
numerical results obtained from CC method with a time step 0.04. The dotted line is the power 
expansion of Eq.(6) with the first four terms. 
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